#!/bin/gawk -f
#
# convMRTABLE: script modifies Aqua restraint file by substituting PDB atom 
# names for the original atom names, using the conversion table which is given 
# in the "MR" file (de-facto standard in Protein Data Bank).
# 
# call:    convMRTABLE mrfile='x' [chain=num] [fixorig=0/1] restraint_file
# example: convMRTABLE mrfile=prot.mr prot.noe
#          convMRTABLE mrfile=dimer.mr chain=90 dimer.noe
# summary message on file: logfil
# output on file: outfil
#
# NOTE: currently only the atom name is converted using the mapping table;
# residue name and number are copied from the restraint (possibly with the
# use of 'chain' correction, see below).
#
# If 'chain' is non-zero, it is assumed that multiple chains occur.
# Residue numbers greater than 'num' (the value of 'chain') will be
# shifted, i.e. the script assumes that the first chain has residue numbers
# up to AND INCLUDING 'num', and the second chain starts again at residue
# number 1.
# The chain identifiers are assumed to be A and B.
# Note: this approach is fairly simple. It will fail for trimers.
# But the whole concept of this script is flawed anyway....
#
# CHAIN directives in the MR table are now also accepted. These should be
# of the form: 
# CHAIN A:
# i.e. a single line appearing before each chain. The chain names given in
# these directives are used to extend the table name definitions, so a map
# on restraints which include chain identifiers can be made.
# The chain variable (if not zero) will now (also) be used to remap residue
# numbers in the MR table. If chain > 0 the mechanism is the same as
# described above. If chain < 0, it is assumed that the residues are numbered
# consecutively, and the effective value of chain is the last residue number
# read before encountering the CHAIN directive.
#
# If 'fixorig' is true, the "orig" atom names are changed from "prefix"
# notation ( 1HB ) to "postfix" notation ( HB1 ).
# This is an option to correct for an inconsistent MR conversion table.
# ===========================================================================
# Jurgen Doreleijers, 30-jan-95
# Ton Rullmann,       12-aug-96; 30-jan-97
# ===========================================================================
# Input and output format: see aqua/doc/restraint_format.txt
# MR file format: see .......
#

BEGIN {
    outfil = "convMRTABLE.tmp"
    logfil = "convMRTABLE.log"

# Store the conversion table
    Restraints_Number = 0
    Notfound_tot = 0
    Notfound_ps = 0
    Notfound_heavy = 0
    Notfound_amide = 0
    Found = 0
    EndModel = 0
    Atom_Number = 0
    Chain_in_mr = 0
    Chain_name = ""
    Table_consec = 0
    if ( chain < 0 ) {
	Table_consec = 1
	chain = -chain
    }

    if ( mrfile == "" ) {
	print "usage: convMRTABLE mrfile='x' [chain=num] restraint_file"
	print "output on file " outfil
	print "log on file " logfil
	exit (1)
    }
# Does the mrfile exist?
    if ( system("test -r " mrfile) != 0 ) {
	print "* mrfile: " mrfile " not found" > logfil
	exit (1)
    }

    while ( ( getline < mrfile ) && !EndModel ) {
	if ( Found ) {
# Does the MR table contain CHAIN directives?
	    if ( match ( $0, /^ *CHAIN *[A-Z]: *$/ ) ) {
		Chain_name = substr( $2, 1, 1 )
		Chain_in_mr = 1
		if ( chain ) {
##		    print "* MR file has CHAIN information" > logfil
##		    print "* conflict with use of input variable 'chain'" >> logfil
##		    exit( 1 )
		    Chain_in_mr = 2
## standard mode: correct residue numbers in MR table using 'chain' value
		    table_chain = chain
## alternate mode: the residue numbers in the MR table are consecutive
                    if ( Table_consec )
			table_chain = Res_num_read
		}
	    }
	    else if ( NF != 10 ) {
# Maybe the end of this model ?
		if ( NF == 0 ) {
		    EndModel = 1
		    break
		}
	    }
	    if ( NF == 10 ) {
		Atom_Number++
		if ( $1 != Atom_Number ) {
		    print "* Atom number: " Atom_Number " not correct at line:" > logfil
		    print "* \"" $0 "\"" >> logfil
		    exit (1)
		}
# Change prefix-notation of "orig" atom name to postfix-notation
		Atom_name = $5
		Res_num = $7
		Res_num_read = $7
		if ( fixorig ) {
		    if ( match( Atom_name, /^[123]/ ) ) {
			seq = substr( Atom_name, 1, 1 )
			sub( /^[123]/, "", Atom_name )
			Atom_name = Atom_name seq
		    }
		}
		if ( Chain_in_mr == 2 )
		{
		    if ( Res_num > table_chain ) Res_num -= table_chain
		}
		Name_atom_residue = Atom_name "_" Res_num "_" Chain_name
# store the PDB atom name, residue name and residue number in three tables
		MRTABLE[Name_atom_residue] = $2
		MRTABLE2[Name_atom_residue] = $3
		MRTABLE3[Name_atom_residue] = $4
	    }
	}
	if ( match( $0, /Start of MODEL/ ) || match( $0, /START_OF_MODEL/ ) )
	    Found = 1
    }
    if ( !Atom_Number ) {
	print "* No atoms found in conversion table" > logfil
	exit (1)
    } 
    print "* " Atom_Number " atoms found in conversion table" > logfil
}

function Map( Atom_name, Residue_nam, Residue_num, Chain_name ) {
# Retrieve the PDB name and residue number for this Name_atom_residue
    Name_atom_residue =  Atom_name "_" Residue_num "_" Chain_name
    if ( MRTABLE[Name_atom_residue] == "" ) {
# Warn if the name is not found, except for pseudo atoms, and keep original name
	if ( !match( Atom_name, /[#*]$/ ) && !match( Atom_name, /^Q/ ) ) {
	    print "* Name of atom and residue not present in table: " Name_atom_residue
	    if ( !match( Atom_name, /^H/ ) )
		Notfound_heavy++
	    if ( Atom_name == "HN" || Atom_name == "H" )
		Notfound_amide++
	}
	else {
	    Notfound_ps++
	}
	Notfound_tot++
	ATM = Atom_name
	RNA = Residue_nam
	RNU = Residue_num
    } 
    else {
	ATM = MRTABLE[Name_atom_residue]
	RNA = MRTABLE2[Name_atom_residue]
	RNU = MRTABLE3[Name_atom_residue]
    }
}
    
# Now do the actual work

NR == 1 {
    print $0 > outfil
    next
}

NR == 2 {
    Count = $2
    print $0 " mr" >> outfil
    next
}

$1 ~ /ULD$/ || $1 ~ /UPLO$/ || $1 ~ /UPP$/ || $1 ~ /LOW$/   {
    Restraints_Number++
    if ( noconv ) {
	print $0 >> outfil
	next
    }
    if ( $0 ~ /INSERT/ ) {
	print "* ERROR: The restraint file has INSERT identifiers." >> logfil
	print "* The conversion table of the MR format does not handle these very well." >> logfil
	print "* From now on the name conversion will be skipped! Sorry..." >> logfil
	noconv = 1
	print $0 >> outfil
	next
    }
    if ( $0 ~ /CHAIN/ ) {
	if ( !Chain_in_mr ) {
	    print "* ERROR: The restraint file has CHAIN identifiers." >> logfil
	    print "* The conversion table of the MR format does not handle these very well." >> logfil
	    print "* From now on the name conversion will be skipped! Sorry..." >> logfil
	    noconv = 1
	    print $0 >> outfil
	    next
	}
	printf( "%-7s ", $1 ) >> outfil
	Map( $6, $4, $5, $3 )
	printf( "CHAIN %s %-4s %-3s %-5s ", $3, $4, $5, ATM ) >> outfil
	Map( $11, $9, $10, $8 )
	printf( "CHAIN %s %-4s %-3s %-5s ", $8, $9, $10, ATM ) >> outfil
	if ( NF == 12 ) 
	    printf( "%10.3f\n", $12 ) >> outfil
	else if ( NF == 13 ) 
	    printf( "%10.3f %10.3f\n", $12, $13 ) >> outfil
	else 
	    printf( "%10.3f %10.3f %10.3f\n", $12, $13, $14 ) >> outfil
    }
    else {
	printf( "%-7s ", $1 ) >> outfil
	if ( !chain ) {
	    Map( $4, $2, $3, "" )
#	    printf( "%-4s %3i %-5s ", RNA, RNU, ATM ) >> outfil
#	    printf( "%-4s %3i %-5s ", $2, RNU, ATM ) >> outfil
#	    printf( "%-4s %3i %-5s ", $2, $3, ATM ) >> outfil
	    printf( "%-4s %-3s %-5s ", $2, $3, ATM ) >> outfil
	}
	else {
	    CHN = "A"
	    if ( $3 > chain ) {
		$3 -= chain
		CHN = "B"
	    }
	    if ( Chain_in_mr )
		Map( $4, $2, $3, CHN )
	    else
		Map( $4, $2, $3, "" )
#	    printf( "CHAIN %c %-4s %3i %-5s ", CHN, $2, $3, ATM ) >> outfil
	    printf( "CHAIN %c %-4s %-3s %-5s ", CHN, $2, $3, ATM ) >> outfil
	}
	if ( !chain ) {
	    Map( $7, $5, $6, "" )
#	    printf( "%-4s %3i %-5s ", RNA, RNU, ATM ) >> outfil
#	    printf( "%-4s %3i %-5s ", $5, RNU, ATM ) >> outfil
#	    printf( "%-4s %3i %-5s ", $5, $6, ATM ) >> outfil
	    printf( "%-4s %-3s %-5s ", $5, $6, ATM ) >> outfil
	}
	else {
	    CHN = "A"
	    if ( $6 > chain ) {
		$6 -= chain
		CHN = "B"
	    }
	    if ( Chain_in_mr )
		Map( $7, $5, $6, CHN )
	    else
		Map( $7, $5, $6, "" )
#	    printf( "CHAIN %c %-4s %3i %-5s ", CHN, $5, $6, ATM ) >> outfil
	    printf( "CHAIN %c %-4s %-3s %-5s ", CHN, $5, $6, ATM ) >> outfil
	}
	if ( NF == 8 ) 
	    printf( "%10.3f\n", $8 ) >> outfil
	else if ( NF == 9 ) 
	    printf( "%10.3f %10.3f\n", $8, $9 ) >> outfil
	else 
	    printf( "%10.3f %10.3f %10.3f\n", $8, $9, $10 ) >> outfil
    }
    next
}

{
    print $0 >> outfil
}

END {
    print "* Processed: "  Restraints_Number " restraints from " FILENAME >> logfil
    if ( Count != Restraints_Number )
	print "* Does not match count of " Count " in file header!" >> logfil
    print "* Not found in table: " Notfound_tot " atoms (used original atom name instead)" >> logfil
    print "*            of which " Notfound_ps " pseudo-atoms, " Notfound_heavy " non-hydrogens, " Notfound_amide " amide protons">> logfil
    print "*  -> potential problem for " (Notfound_tot-Notfound_ps-Notfound_heavy-Notfound_amide) " protons" >> logfil
}
